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Abstract 

We propose a method for post-processing an ensemble of multivariate forecasts in order to 
obtain a joint predictive distribution of weather. Our method utilizes existing univariate post- 
processing techniques, in this case ensemble Bayesian model averaging (BMA), to obtain es- 
timated marginal distributions. However, implementing these methods individually offers no 
information regarding the joint distribution. To correct this, we propose the use of a Gaussian 
copula, which offers a simple procedure for recovering the dependence that is lost in the esti- 
mation of the ensemble BMA marginals. Our method is applied to 48-h forecasts of a set of 
five weather quantities using the 8-member University of Washington mesoscale ensemble. We 
show that our method recovers many well-understood dependencies between weather quanti- 
ties and subsequently improves calibration and sharpness over both the raw ensemble and a 
method which does not incorporate joint distributional information. 



1 Introduction 



Mesoscale weather forecasting is usually conducted via a forecast ensemble where the ensemble 
members differ by the boundary conditions and/or the parameterization of the model physics in 
the numerical weather prediction (NWP) model (Leutbecher and Palmer] 2008| ). While ensemble 
systems aim to reflect and quantify sources of uncertainty in the forecast, they tend to be biased, 



and they are typically underdispersed ( |Hamill and Colucci| |1997| ). The predictive performance of 
an ensemble forecast can thus be substantially improved by applying statistical postprocessing to 
the model output in which the forecasts are corrected in coherence with recently observed forecast 
errors (Wilk sand Hamill| |2007| ). It has further been argued that the statistical postprocessing tech- 
niques should be probabilistic in nature and return full predictive distributions, see e.g. IGneiting 



and Raftery] ( |2005[ ). State of the art approaches of this type include ensemble model output statis- 



tics or nonhomogeneous Gaussian regression (Gneiting et al. , 2005 ; Thorarinsdottir and Gneiting 
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2010| |Thorarinsdottir and Johnson], |2011|) and kernel dressing or ensemble Bayesian model aver- 



aging dRaftery et alii |2005| [Fortin et al.j [20061 |Sloughter et aL| |2007l |20T0l ). 

However, many of these techniques consider only a single weather quantity at a fixed look- 
ahead time, without taking potential spatial dependencies into account. The postprocessed forecasts 
may thus violate the multivariate correlation structure of the original ensemble forecasts and the 
observations. For low-dimensional multivariate settings, the correlation structure can be modelled 
directly with a parametric model. Parametric postprocessing techniques that model the spatial 
dependencies explicitely through a geostatistical model have e.g. been proposed for temperature 
(Berroca l et al.[ |2007) and precipitation ( |Berrocal et al.H200"8j ). Similarly, [Pinson ( 2011| ), [Slou ghter 



et al.| ( [20TT] ), and |Schuhen et al.| ( |2012[ ) consider bivariate probabilistic forecasts for wind vectors. 



In higher dimensions, joint parametric modelling becomes cumbersome, especially when the 
marginal distributions are assumed to be of different types as is e.g. the case for temperature, 
precipitation, and wind speed. For such data, copula methods are advantageous as they allow for 
independent modelling of the marginal distributions and the multivariate dependence structure of 
the rank statistics, see e.g. Genest and Favre ( |2007 ). We propose a multivariate postprocessing 
framework where in a first step, established Bayesian model averaging (BMA) methodology is 
applied independently to each weather variable to obtain calibrated and sharp marginal predictive 
distributions. In a second step, the marginal distributions are conjoined in a multivariate framework 
using a Gaussian copula model, see e.g. Hoff|(2007). 

We apply our method to 48-hour ahead forecasts of five weather variables obtained from the 
eight-member University of Washington mesoscale ensemble ( |Eckel and Mass] |20Q5j ) at 60 ob- 
servation locations in the North American Pacific Northwest in 2008. The variables we consider 
are daily maximum and minimum temperature, sea level pressure, precipitation accumulation, and 
maximum wind speed. The marginal predictive distributions of temperature and pressure are de- 
fined on the entire real axis while wind speed takes values on the positive real axis only. Fur- 
thermore, the marginal predictive distribution for precipitation accumulation takes values on the 
non-negative real axis with an additional point mass in zero. While estimation of the marginal 
distribution requires incorporation of these features, the Gaussian copula remains largely agnostic 
to such factors, highlighting its flexibility when working distributions with mixed marginals (Hoff , 
2007[ ). An example of the resulting multivariate predictive distribution is given in Figure [TJ 



Copula methods are widely used for prediction problems in hydrology, see e.g. |Genest and 



Favre| ( |2007l ), |Scholzel and Fr iederichs (2008 ), [Kao and Govindaraju| ( |2010| ) and references therein. 



In the context of multivariate statistical postprocessing, two discrete copula approaches have been 
proposed where the multivariate rank structure is inherited from either past observations or the 
ensemble forecasts. The Schaake Shuffle (Clark et al. 2004[ ) reorders a postprocessed ensemble 
based on past observations in order to recover the multivariate rank structure in the observed data. 
Ensemble copula coupling (ECC) ( |Schefzik||2011[ ), on the other hand, constructs multivariate sam- 
ples from the marginal predictive distributions in such a way that the multivariate rank structure 
of the original NWP ensemble is preserved. Our work, by contrast, uses a Gaussian model for 
dependence, which requires estimation of a single model component (the multivariate correlation 
matrix). We feel this approach has the potential to scale considerably better than approaches based 
on estimation of multivariate ranks, which would require a considerable number of observations in 
high dimensions to estimate properly. Note that the method proposed by Pinson (2011 1 may also 
be considered a copula approach. 

The paper is organized as follows. Description of the data and the methods are given in Sec- 
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tion|2j We shortly review the univariate BMA approaches and the multivariate verification methods 
that we apply and give a detailed description of the Gaussian copula approach. The results of the 
case study are presented in Section [3] and the paper ends with conclusions in Section |4} 



2 Data and methods 
2.1 Data 

In our case study, we employ daily 48-h forecasts based on the University of Washington mesoscale 
ensemble (UWME; Eckel and Mass 2005) with valid dates in the calendar year 2008. The UWME 
is an eight-member multianalysis ensemble which then was based on the Fifth-Generation Penn 
State/NCAR Mesoscale Model (MM5) with initial and lateral boundary conditions obtained from 
operational centers around the world. Currently, the UWME uses the WRF mesoscale model. 
Further information as well as real time forecasts and observations can be found on the website 
|http : / /www, atmos .Washington . edu/~ens/uwme . cgij 

The forecasts are made on a 12 km grid over the Pacific Northwest region of Western North 
America. To obtain a forecast at a given observation location, the forecasts at the four surrounding 
grid points are bilineraly interpolated to that location. We consider observation locations in the US 
states of Washington, Oregon, Idaho, California, and Nevada, see Figure [3] The daily observations 



are provided by weather observation stations in the Automated Surface Observing Network (Na 



tional Weather Service} |1998| ). We consider 2-m maximum and minimum temperature, sea level 



pressure, 10-m maximum wind speed, and 24-h precipitation accumulation. Forecasts and obser- 
vations are initialized at 00 UTC which is 5pm local time when daylight saving time operates and 



4pm local time otherwise. Quality control procedures as described by Baars (2005) were applied 
to the entire data set, removing dates and locations with any missing forecasts or observations. 

For the calendar year 2008, we consider 60 distinct observation locations which have between 
95 and 271 days in which all ensemble forecasts and verifying observations were available. Addi- 
tional data from 2006 and 2007 were used to provide an appropriate rolling training period for all 
days in 2008 and to estimate the multivariate correlation structure. 

2.2 Ensemble Bayesian model averaging 

Bayesian model averaging (BMA) was originally developed as a method to combine predictions 
and inferences from multiple statistical models ( Learner} |1 97 8| ). Raftery et al. ( 2005| ) extended the 



use of BMA to statistical postprocessing for forecast ensembles. In this context, the method is 
a kernel dressing approach where each ensemble member Xk is associated with a kernel density 
9k(y\xk)- The ensemble BMA predictive density is then given by a mixture of the individual kernel 
densities, 



A 



f(y\x 1} ...,x K ) = ^2<ij k g h (y\x k ), (1) 



k=l 



where the weights uo k are assumed to be non-negative with Y^=i u k = 1- The choice of the kernel 



g k depends heavily on the weather variable of interest: Raft ery et al.| ( |2005[ ) consider temperature 
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Table 1 : The ensemble BMA kernel functions for the different weather variables and the associated 
link functions for each mean value and variance. The gamma distribution is parameterized in terms 
of shape and scale with mean a/3 and variance a/3 2 . 



Variable 


Range 


Kernel 


Mean 


Variance 


Temperature 


y e R 




bok + b± k x k 


a 2 


Pressure 


y e R 




bok + b\ k x k 


a 2 


Wind speed 


yeR + 


T(a k ,(3 k ) 


bok + b\ k x k 


c + c x x k 


Precipitation amount 


y 1 /' 3 e R- 


h T(a k ,/3 k ) 


bok + b lk xf 


Co + CiX k 



and pressure for which Gaussian kernels seem appropriate, while Sloughter et al. ( 2010| ) apply 
gamma kernels to wind speed forecasts. 

Precipitation has to be treated slightly differently as precipitation observations are non-negative 



with a large number of zero observations. Sloughter et al. (2007) propose a solution to this where 
the kernel density g k is modelled in two parts, 

g k {y\x k ) = F{y = 0\x k )l{y = 0} + F{y > 0|a?k)/i*(2/ 1/3 |a; fc )l{2/ > 0}, 



where h k is a gamma density and 



P(y = o|/ fc 



exp(q ofc + ai k xl /s + a 2 kh) 
aikxl /3 + a 2 kh) 



exp(a fc 



with 5k = 1 if x k = and 5 k = otherwise. Note that h k is a predictive density for the cube root 
of the precipitation amount. However, the resulting probabilistic forecast can easily be expressed 
in terms of the original amounts. 

Table[T]gives an overview over the different BMA models we consider, as well as the associated 
link functions for the mean value and the variance of each kernel density. Further variants of the 



ensemble BMA method which will not be considered here include the work by Roquelaure and 



Bergot (2008 1, Bao et al. ( |2010 ), and[Chmielecki and Raftery (2010). For the parameter estimation, 
we apply the R package ensembleBMA which provides estimation methods for all the models 
listed in Table [T|( R Development Core Team[|20imFraley et al.[[2011| ). 



2.3 Gaussian Copulas 



The ensemble BMA methods discussed in Section 2.2 have proven to work well for post-processing 
ensemble forecasts of univariate quantities. However, by using these methods for each quantity 
individually, no attention has been paid to the joint distribution of weather quantities. In this section, 
we outline a Gaussian copula approach that allows us to recover the dependence between weather 
quantities and construct a post-processed joint distribution. 

Suppose that we have p weather quantities of interest, with marginal distributions Fi, . . . , F p , 
where 



Fj(y) 



fj{u\xij, . . .,x K j)du. 



(2) 
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In Equation [5J fj(u\xij, . . . , XKj) represents the ensemble BMA density discussed in Equation [T] 
for variable j, evaluated at u and depending on the ensemble members x±j, . . . , XKj- Now let C 
be a p x p correlation matrix, i.e. a positive definite matrix with unit diagonal. Under a Gaussian 
copula, the joint distribution F of the weather quantities takes the following form 

F( Vl , . . .,y p \C) = $ P ($- 1 (F 1 ( Z/1 )), • • • , ®-\F p (y p ))\C), (3) 

where is the inverse cdf of a standard Gaussian distribution and is the cdf of a 

p— variate Gaussian distribution with covariance matrix £. The Gaussian copula is a particularly 
tractable type of copula model, as it requires only the marginal distributions Fi, . . . ,F P and the 



correlation matrix C to be fully defined (see |Nelsen| ( |2006[ ) for a more detailed account of general 
copula models). 

The Gaussian copula lends itself to a useful construction. Let 

Z ~ JV p (0, C) 



and for j = 1 , . . . , p set 
where 



Yj = F~ 1 (^(Zj)), 



F j l (u) = sup{y : F d (y) < u} 

denotes the psuedo-inverse of the marginal Fj. Then we have that Y = (Yi, . . . , Y p ) ~ F. The 
construction also highlights that each Yj is marginally distributed according to Fj. 

The construction above thus creates a link between an observation Y sampled from F and a la- 
tent Gaussian factor Z. In particular if Fj is a continuous marginal distribution we can immediately 
see that 

Zj = ^(FjiYj)). 

This indicates that for the majority of weather quantities, which have fully continuous distributions, 
given Fj and an observed yj we can directly infer a latent Zj. In the case of precipitation, the 
situation is slightly more nuanced. In general, suppose that Yj G [0, +00) where Fj(0) = a with 
< a < 1 and Fj is otherwise continuously increasing on (0, +00). Then we have that 

-00 < Zj < <S>-\a) 

Zj = 

when Yj > 0. If we therefore collect several observations y^, . . . , y( T ) we may infer latent Gaus- 
sian observations zW, . . . , and thereby estimate the matrix C. 

Our process for forming the ensemble post-processed joint distribution of weather quantities 
builds on the logic above. Suppose that we have a collection . . . , y( T ) of observations over T 
days. We assume that each 

y M~F (t) HC) 



when Yj = and 



where 



F« („?>, . . . ,y(«)|c) = $ P (ff(y?>)) , • • • , ^ {F}%®)) |C 
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and F^ denotes the ensemble BMA marginal distribution for weather quantity j at timepoint 

t using the K ensemble members xf\ . . . x^ at timepoint t. This framework associates each 
observation y(*' with its own Gaussian copula F^, but all copulas share one correlation matrix C. 
The assumption of a common correlation matrix C can be relaxed, as we outline in the Section |4} 

Using y(*) and the marginal distributions F-J , . . . , Fp we may then infer a latent Gaussian 
as discussed above. 

While each was given a separate distribution we note that the associated zW are all 
distributed J\f p (0, C). We therefore use these latent Gaussian observations to estimate C, for in- 



stance by taking the sample correlation matrix. |Hoff (2007) discusses more involved methods for 



estimating C, we discuss the inclusion of these methods in Section HI 

Now consider forming a predictive distribution for the timepoint s (coming sometime after T), 
based on the K ensemble members x^ , . . . , xj^ and the estimate C. Our method proceeds by first 

(s) (s) 

forming the ensemble BMA predictive marginals F± , . . . , Fp ' and then setting 

Y (s) ~ F (s) (-|C). 

While this joint predictive distribution may not have an easy analytic structure, a sample Y may be 
obtained by first sampling 

Z ~ JV p (0, C) 

and then setting 

By sampling a large number of Y in this manner, we are able to effectively describe the entire joint 
predictive distribution. As noted above, the marginal distributions for each individual quantity of 

(si (s) 

this sample remain the ensemble BMA marginals F± , . . . , Fp . 



2.4 Multivariate forecast verification 



To assess the quality of the multivariate forecasts, we apply the methods described in Gneiting 
et al. ( 2008[ ). For inspecting calibration, we use the multivariate rank histogram (MRH) which is a 



direct generalization of the univariate verification rank histogram or Talagrand diagram (Anderson 



1996; [Hamill and Colucci 1997, Talagrand et al. , 1997). The only challenge lies in defining a 
multivariate rank order, as no natural ordering exists for multivariate vectors. Here, we use the 
multivariate ordering described in Gneiting et al.| (2008). 

A forecast is said to be calibrated if the resulting MRH is close to being uniform. To quantify 
the deviation from uniformity, we use the discrepancy or reliability index A, 



A 



m+l 



m + l 



where Q is the observed relative frequency of rank j (Delle Monache et al. 2006) 



The sharpness of a univariate predictive distribution or an ensemble forecast can easily be as- 
sessed by the corresponding standard deviation. In the multivariate case, we employ a generaliza- 
tion of the standard deviation, the determinant sharpness (DS), 



DS = (det S 
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where S is the covariance matrix of an ensemble or a multivariate predictive distribution for a d- 
dimensional quantity. For ensemble forecasts, the matrix is generated by empirical variances and 
correlation of the ensemble. 

Scoring rules for the verification of deterministic or probabilistic forecasts are well known and 
have been widely used in forecast assessment. We consider multivariate extensions of the absolute 



error and the continuous ranked probability score (Matheson and Winkler 1976 Hersbach, 2000). 



The absolute error generalizes to the Euclidean error (EE), 

EE(F,y) = ||A*-y||, 

where n is the median of F. For an ensemble or a sample from a continuous distribution, fi is 
defined as the vector that minimizes the sum of the Euclidean distance to the individual forecast 
vectors, 

mm 



The vector n can be determined numerically using the algorithm described in Vardi and Zhang 
( |2000p as implemented in the R package ICSNP 



For a generalization of the continuous ranked probability score, |Gneiting and Raftery] ( |2007| ) 
introduce the energy score (ES), 

ES(F,y)= J B J r||X-y||-^ F ||X-X , ||, 

where 1 1 • 1 1 denotes the Euclidean norm and X and X' are independent random vectors with distri- 
bution F. If F is the cumulative distribution function associated with a forecast ensemble of size 
m, the energy score can be computed as 



1 rn 1 m rn 

ES Wy) = ^Eii^--yii-wEE 



2m 2 

3=1 t=i 

Generally, the energy score may be approximated by 

- n 1 ™ 

ES(F,y)«- £ || Xi -y||-- £ 



/ 1 



where {xj}" =1 and {x^}™ =1 are two independent samples from F. 

We assign the forecasting methods a score by averaging the scoring rules over all locations and 
time points in the test set. Both the energy score and the Euclidean error are negatively oriented 
such that a smaller score indicates a better predictive performance. The values of the weather 
variables we consider are given on scales which vary by several orders of magnitude. For this 
reason, we normalize the components before we calculate the scores using the observed mean 
values and standard deviations over the test set. 



3 Results 

We now give the results of applying the established BMA framework together with our multivariate 
Gaussian copula approach to maximum and minimum temperature, sea level pressure, maximum 
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Table 2: Estimated correlation matrix at the station KSEA, Sea-Tac Airport, based on data from the 
calendar year 2007. 





maxwsp 


precip 


mintemp 


maxtemp 


pressure 


maxwsp 


1 


-0.016 


0.032 


0.139 


-0.123 


precip 


-0.016 


1 


-0.001 


-0.174 


-0.015 


mintemp 


0.032 


-0.001 


1 


0.239 


-0.110 


maxtemp 


0.139 


-0.174 


0.239 


1 


-0.203 


pressure 


-0.123 


-0.015 


-0.110 


-0.203 


1 



wind speed, and precipitation over the North American Pacific Northwest in 2008. See Section 2.1 
for a detailed description of the data. The BMA univariate postprocessing is applied at each obser- 
vation location separately. Based on an exploratory analysis using a subset of the data set, we use a 
40-day sliding training period for the parameter estimation. That is, the training period consists of 
the 40 most recent days prior to the forecast for which ensemble output and verifying observations 
were available. In terms of calendar days, this period typically corresponds to more than 40 days. 

3.1 Results at Sea-Tac Airport 

To show the behavior of our methodology in-depth, we first focus on the KSEA observation sta- 
tion, located at Sea-Tac Airport, a major transportation hub in the area. Using all available data 
from 2007, we run the ensembleBMA methodology for each of the five variables as described in 



Section 2.2 We then use the observations for these data and the estimated marginal distributions to 



infer a latent vector z as described in Section 2.3 This is performed separately for each day in 2007 



and the resulting latent data are then used to estimate a single correlation matrix. Table [2] shows the 
entries of this correlation matrix. 

The correlations estimated in Table[2]show several clear patterns of interaction. We see a strong 
negative correlation between pressure and both maximum and minimum temperatures. This is in 
line with the understood inverse relationship between temperature and pressure systems. We also 
see an intuitive positive correlation between the minimum and maximum temperatures. 

This estimated correlation matrix is then carried forward into an analysis of 2008, which we use 
for verification. For each day with observations in 2008, we again run ensembleBMA individually 
for each weather quantity to obtain estimated marginal predictive distributions. We then use the 
estimated correlation matrix and these marginal distributions to obtain 20,000 samples from the 



joint predictive distribution for each day, as described in Section 2.3 

Figure [T] shows a pairwise plot of this joint predictive distribution, for the date January 1st, 
2008. For each pair of variables, the figure shows a heatmap indicating regions of significant 
probability mass - with lighter regions being higher values - as well as points showing the 8 
ensemble members (circles) and the observed level (a square). The marginal predictive distribution, 
given by the ensembleBMA methodology, is shown for each variable along the diagonal. 

In this figure we can see that the correlation structure in Table [2] has been carried over to the 
predictive distribution: the positive correlation between maximum and minimum temperatures is 
evident as well as the negative correlation of each of these quantities with pressure. Further, we 
can see the effect of post-processing, as the predictive distributions are often centered away from 
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Figure 1: Estimated joint predictive distribution for January 1st, 2008 at the KSEA observation 
station along with ensemble predictions (circles) and verifying observation (square). In the pair- 
wise plots, lighter areas correspond to regions of higher probability mass. The diagonal shows 
the marginal predictive distribution for each quantity. Wind speed is given in meters per second, 
precipitation in millimeters, temperature in degrees C, and pressure in millibars. 

2468 12 01234 579 11 68 12 16 1006 1010 1014 




468 12 01234 579 11 68 12 16 1006 1010 1014 



the ensemble members (indicating the bias correction properties) and display greater spread in the 
distribution than is evident in the ensemble. 

Table [3] shows that according to a number of verification metrics, the Gaussian copula ap- 
proach yields an improvement in predictive performance. The table compares the Gaussian copula 
approach to an alternative where dependence is not modeled, which we call the independence ap- 
proach, as well as the raw UWME ensemble averaged over 2008 for the KSEA observation station. 
Note that the Gaussian copula and independence approaches have the same marginal distributions, 
and thus differ only in the manner in which the joint distribution is constructed. 

Both the copula and independence approaches improve considerably on the raw ensemble in all 
metrics except the determinant sharpness (DS). This is not surprising since the 8 members of the 
raw ensemble will impart greater sharpness at the expense of calibration. We also see that the copula 
approach improves on the independence approach for all metrics. The values of the reliability index 
and the DS show that the predictive distribution for the copula approach is both better calibrated 
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Table 3: Predictive performance of the copula methodology, the independence approach, and the 
raw UWME ensemble at the KSEA observation station. Results are averaged over the 271 days in 
2008 for which forecasts and verifying observations were available. The performance is measured 
by the energy score (ES), the Euclidean error (EE), the reliability index (A), and the determinant 
sharpness (DS). 





ES 


EE 


A 


DS 


UWME 


0.938 


1.081 


0.185 


0.566 


Independence 


0.637 


0.982 


0.047 


7.516 


Copula 


0.636 


0.982 


0.019 


6.971 



(a) UWME (b) Independence (c) Copula 

Figure 2: Multivariate rank histograms for the copula and independence approaches as well as the 
UWME ensemble for the KSEA observation station over the 295 days in 2008 for which forecasts 
and verifying observations were available. 



and somewhat sharper than for the independence approach, a consequence of using a non-diagonal 
correlation matrix. The Euclidean scores are essentially the same for the two approaches, since 
they can largely be expected to return similar median values. The combination of similar median 
and improved sharpness leads the copula approach to have a lower energy score. 

Figure [2] shows the multivariate rank histograms for the two approaches, as well as for the 
raw ensemble. The figure shows that both methods improve calibration considerably over the 
raw ensemble. However, as shown in the figure, in the independence approach the final bins are 
somewhat less filled than under the copula approach, though neither returns a completely uniform 
rank histogram. 

3.2 Results over the Northwest US 

A similar analysis as that above was run for 60 separate observation stations in the Northwest US. 
First, ensembleBMA was run individually for each station, day and weather quantity during the 
period of 2007. Verifying observations were then used to estimate a correlation matrix separately 
for each observation station. While these estimates were performed locally, there is considerable 
agreement in estimated correlations between individual stations. Figure [3|a) shows the pairwise 
estimated correlation between minimum and maximum temperature at each observation station. 



10 



This plot shows that the majority of estimates are positive, as expected, with correlations as high as 
0.51. 

Figure ^b) shows a similar plot, but with the correlation between minimum temperature and 
pressure plotted for each observation station. While the previous figure showed some similarity in 
estimated correlations, this figure shows a considerable agreement between observation stations. 
We see that all values are quite negative, with the majority of estimates between —0.2 and —0.4. 
Furthermore, there appears to be a spatial component to the estimates. The values near the top of 
the Puget Sound are all roughly between —0.1 and —0.2, while those closer to the Seattle/Tacoma 
area are grouped between —0.22 and —0.26. A tight group of observation stations in the Columbia 
River Valley, on the border of Washington and Oregon, all have correlations between —0.31 and 
—0.35 and finally those in Eastern Washington and Eastern Oregon exhibit stronger correlations 
typically below —0.4. Figures [3fa) and ^b) therefore suggest that an important feature of the joint 
distribution is captured through the Gaussian copula methodology that is ignored in the indepen- 
dence approach. 
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(b) Minimum temperature and pressure 



Figure 3: Estimated correlation between (a) minimum and maximum temperature, and (b) mini- 
mum temperature and pressure at 60 observation stations in the Northwest US using 2007 data and 
the Gaussian copula approach. 

Figure[4]compares the multivariate rank histograms (averaged over all observation stations) for 
the copula and independence approaches as well as the raw UWME ensemble. We see a similar 
result across all stations as for the KSEA station above: namely that the highest bin is under- 
occupied in the independence approach in comparison to the copula approach. Furthermore, both 
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(a) UWME 



(b) Independence 



(c) Copula 



Figure 4: Multivariate Rank Histograms for the copula and independence approaches as well as 
the raw UWME ensemble taken over all available observations at 60 observation stations in the 
Northwest US over the year 2008. 



Table 4: Predictive performance of the copula and independence approaches as well as the raw 
UWME ensemble. Results are averaged over 60 observation stations in the Northwest US and 
all days in 2008 for which forecasts and verifying observations were available. The predictive 
performance is measured by the energy score (ES), the Euclidean error (EE), the reliability index 
(A), and the determinant sharpness (DS). 





ES 


EE 


A 


DS 


UWME 


0.881 


1.061 


0.161 


0.811 


Independence 


0.586 


0.914 


0.071 


1.945 


Copula 


0.585 


0.914 


0.066 


1.905 



methods improve the calibration considerably compared to the raw ensemble. 

Table [4] presents verification scores averaged over all the 60 stations and all available days in 
2008. The results here are broadly consistent with those reported for KSEA. We see that the de- 
terminant sharpness is improved in the copula versus independence approach, while the Euclidean 
score is essentially the same. These two factors lead to an improvement in the energy score in the 
copula approach. Calibration also appears to be improved, as shown through lower energy score 
and reliability index for the copula approach. 



4 Conclusions 

We have proposed a method for constructing a joint predictive distribution based on post-processing 
an ensemble of forecasts for multiple weather quantities. Our method utilizes existing techniques 
for post-processing univariate weather quantities and then leverages a Gaussian copula to tie these 
individual marginal distributions together. The method is relatively simple, it requires little ad- 
ditional computational effort after the univariate marginals are formed and it retains the marginal 
distributions learned from ensemble BMA. We have then shown that the method yields a calibrated 
and sharp distribution, using data from the Pacific Northwest. 
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In this paper, we focused on using the ensemble BMA methodology to construct marginal pre- 
dictive distributions. This was done since the methodology has already been established for several 
interesting weather quantities and software exists to implement these methods. However, in prac- 
tice any method for forming marginal distributions could replace the ensemble BMA framework 
and essentially nothing of the overall methodology would be affected. This highlights the flexibility 
of the copula approach. 

Alternative approaches to constructing joint distributions such as the Schaake shuffle and the 
ECC approach use a discrete copula to construct the joint distribution. These methods learn a mul- 
tivariate rank structure by investigating ranks in historical observations (Schaake) or the ensemble 
(ECC). We feel our approach offers a useful alternative to these methods for two reasons. First, 
the joint distribution is modeled using a Gaussian distribution. This does not require subsequent 
samples from the predictive distribution to strictly obey any observed rank structure and, in our 
opinion, may scale better to high dimensions. The statistics literature is currently investigating 
these questions (see |Hoff et aL ( 201 1| )) and subsequent work should be undertaken to compare 
these methods in the ensemble post-processing context. The second feature of our method is that 
it aims to model the joint distribution after post-processing. By running the ensemble BMA ap- 
proach and then learning the latent Gaussian factor using a verifying observation, we are implicitly 
modeling the joint residual structure implied by the ensemble BMA method itself. This approach, 
in our opinion, more appropriately reflects the modeling process that is being performed. 

Section [23] outlined the basics of Gaussian copula theory and provided a simple algorithm for 
calculating latent Gaussian factors and the subsequent correlation matrix. We presented this ap- 
proach as it is relatively straightforward and captures the main features of interest. More involved 
estimation strategies exist, for instance the Bayesian approach of Hoff (2007 1 and subsequent work 
should consider these frameworks to assess whether they offer improvements in predictive perfor- 
mance. 



As noted in Section 12.31 we assume a constant correlation factor for all observations when es- 
timating the Gaussian copula. Certainly, such an assumption may be an over-simplification. More 
complicated methods that include time-varying correlation factors are possible, but these methods 
are both more difficult to describe and to estimate. Further, preliminary investigations on our part 
revealed essentially no added benefit to considering time-varying correlation models in the exam- 
ple discussed in Section [3} However, it may be useful to reconsider the estimation strategy in the 
future, especially in higher-dimensional situations. 

The Gaussian copula has received occasional criticism (see Mikosch (2006)) for not accurately 
capturing dependence in the tails of multivariate distributions. This criticism is germane, but not 
a central concern to the modeling task undertaken in this paper. Our goal has been to construct a 
multivariate distribution after post-processing an ensemble of forecasts. If extreme weather events 
were possible, it is probable that this information would already be incorporated to some extent 
in the ensemble itself, and the copula would simply indicate the variability about these extreme 
values. Furthermore, the direct modeling task has focused on day-to-day weather prediction, as 
opposed to focusing on extreme events. If the joint modeling of extreme weather is the central 
concern of the forecaster, it is likely that an alternative copula model may be preferable. 

While we have shown promising results in modeling five weather quantities simultaneously, 
the long term goal is undoubtedly to model weather jointly in the spatial domain and for multiple 
variables. Once spatial factors are added, the dimensionality of the model can increase rapidly and 
issues related to spatial covariation estimation are necessary to solve. There have been several re- 



13 



cent advances in using fast computational methods for Gaussian Markov random fields dLmdgren 



et al.[ 2011| ), which could prove useful in constructing high-dimensional joint distributions based 



on post-processed ensemble forecasts. Our current steps in advancing the methods discussed above 
have been to merge this literature with that of ensemble post-processing. 
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